#!/usr/bin/perl -w
# this script was used to help generate test files for CodonFreqStdDev.pl
# it generates a fasta sequence where the distribution of each nucleotides is
# exactly that in $codonTableFile, with some fractional difference $diff
use strict;
use warnings;

if (@ARGV != 2) {
  print "USAGE: $0 <Codon Table> <Fraction of Deviation>\n";
  exit;
}

my ($codonTableFile, $diff) = @ARGV; 

my %codonTable = getCodonTableFrom ($codonTableFile);
my $seq = "";

foreach my $codon (keys %codonTable) {
  if($codon =~ m/^[TC]/) {
    $codonTable{$codon} += $diff*$codonTable{$codon};
  } else {
    $codonTable{$codon} -= $diff*$codonTable{$codon};
  }
  my $count = round($codonTable{$codon} * 10000);
  for (my $i=0; $i<$count; $i++) {
    $seq .= $codon;
  }
}

$seq = lc $seq;
$seq =~ s/(.{1,80})/$1\n/g;
print ">seq1\n$seq\n";

sub round {
  my ($float) = @_;
  return sprintf "%.0f", $float;
}

# parses a codon frequency table from the Codon Usage Database
sub getCodonTableFrom {
  my ($codonTableFile) = @_;
  my %codonTable;
  open (INPUT, $codonTableFile) or die "Can't open $codonTableFile: $!";

  while( my $line = <INPUT> ) {
    if ( $line =~ /^\s*$/ ) {
      next;
    }
    $line =~ s/\(.*?\)//g; # get rid of parens and the stuff inside of them
    $line = uc $line;
    $line =~ tr/U/T/; # convert all U's to T's (we're expecting DNA)
    my @stats = split (/\s+/, $line); # split on whitespace
    for (my $i=0; $i<4; $i++) {
      # divide by 1000 so that the sum of the frequencies is 1
      $codonTable{$stats[2*$i]} = $stats[2*$i+1]/1000; 
    }
  }
  return %codonTable;
}

